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Motivated by recent experimental results, we study beating dark-dark solitons as a prototypical coherent 
structure that emerges in two-component Bose-Einstein condensates. We showcase their connection to dark- 
bright solitons via SO(2) rotation, and infer from it both their intrinsic beating frequency and their frequency of 
oscillation inside a parabolic trap. We identify them as exact periodic orbits in the Manakov limit of equal inter- 
and intra- species nonlinearity strengths with and without the trap and showcase the persistence of such states 
upon weak deviations from this limit. We also consider large deviations from the Manakov limit illustrating that 
this breathing state may be broken apart into dark-antidark soliton states. Finally, we consider the dynamics 
and interactions of two beating dark-dark solitons in the absence and in the presence of the trap, inferring their 
typically repulsive interaction. 



I. INTRODUCTION 



One of the principal themes of study in the emerging field of atomic Bose-Einstein condensates (BECs) is the examination of 
the coherent structures that arise in them (TJ-E). When such explorations started over a decade ago |5 - 9], they were considerably 
hindered by either geometric or thermal effects, which were detrimental towards the lifetime of dark solitons and vortices that 
can be formed in repulsive BECs. Yet, the newer generations of experiments have enabled considerable strides towards the 
observation of dynamics and interactions of such nonlinear waveforms fT0lfT5ll . 

In addition to the above context of single-component BECs, soliton and vortex states may also arise in multi-component 
condensates, such as the two-component pseudo-spinor BECs, or the three and higher component spinor BECs dElE]. One 
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FIG. 1 : (Color online) Prototypical experimental images of dark-bright and dark-dark solitons in a two-component BEC. The two components 
are vertically offset for separate imaging. All dynamics occur with vertically overlapped components before the imaging procedure. Clear 
examples of dark-bright and dark-dark solitons are marked as DB and DD respectively. In the fourth panel, the red (thick) line shows a radially 
integrated cross section of the upper component in the boxed region of the third panel, while the black (thin) line shows the cross section of 
the lower component. The |F, 7tif) hyperfine states used for these images are given to the right of each component. 
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of the prototypical examples of a soliton state in these settings is the so-called dark-bright (DB) soliton (T6J021]. Experimental 
images of DB solitons in a two-component BEC are presented in Fig. [T] The BEC in this figure is comprised of two different 
hyperfine states of 87 Rb, and the solitons are generated by subjecting the BEC to inter-component counterflow; details of this 
technique are described in Refs. [010 [T9ll . In each panel, the atom clouds of the two components are vertically offset for 
imaging only, while all the dynamics leading to the soliton formation occurs in overlapped clouds. Clear examples of dark- 
bright solitons are marked as "DB" in the figure and they consist of a dark soliton in one component that is coupled to a bright 
soliton in the second component. These structures can be thought of as "symbiotic" (or even parasitic) states because their 
bright component cannot be supported alone in the case of repulsive interactions ; in fact, the bright soliton is only sustained 
because of the presence of its dark-counterpart, which operates as an external trapping potential. Although dark-bright solitons 
(and even a prototypical interaction thereof) were first observed some time ago in the context of nonlinear optics l20ll2T1l . their 
observation in recent atomic BECs experiments 1 10 ] triggered a sizeable burst of research activity centered around them. Topics 
of study included (but were not limited to) multi-DB soliton solutions from the viewpoint of integrable systems |22], numerical 
study of DB soliton interactions l23ll , discrete DB solitons l24l . experimental realizations of DB soliton trains |18 ], DB soliton 
oscillations and interactions (25j|26l, as well as interaction of DB solitons with localized impurities l27l . 

Recently, a "cousin" of these DB solitons, namely the dark-dark (DD) soliton — which involves two dark solitons but with 
potentially a breathing oscillation between their densities was also experimentally observed lfT9l . Pertinent examples are marked 
as "DD" in Fig.[T] These solitons show interesting dynamics in which they periodically change their form, from the one shown 
in the first panel to the one shown in the second panel, and back (note the order of the hump/notch features in each of the DD's 
component; see also Fig. [5] below). Such "beating dark-dark solitons" are expected to emerge in the integrable two-component 
(so-called Manakov) limit of the relevant mean-field theoretic models |[28l and were, in fact, earlier observed in numerical 
experiments involving the dragging of defects through the binary condensates (29). 

The current experimental advances, such as the ones leading to the soliton images of Fig. [T] motivate the present theoretical 
study, in which we revisit DD soliton states at the integrable Manakov limit and extract information from their connection to 
the DB solitons (section II). These results are corroborated by the identification of such single DD soliton states, as genuine 
periodic orbits of the Manakov case (with and without a trap) and the study of their stability, internal modes and associated 
near-equilibrium dynamics (section III). In addition, we examine the dynamics of individual such solitons, upon departure from 
the integrable limit (section IV). Experimentally it has also become possible to generate several solitons, and even solitons of 
different types, in a single BEC - see, e.g., the third panel of Fig. [T] which demonstrates the coexistence of dark-bright and 
dark-dark solitons. Although the experimentally exploited counterflow between the two components is beyond the scope of 
our current analysis, these experimental findings motivate our investigation of the interactions between two dark-dark solitons 
(section V). Finally, conclusions of our study, as well as a number of interesting perspectives for future work, are also presented 
(section VI). 



II. DARK-BRIGHT AND DARK-DARK SOLITONS: THEORETICAL BACKGROUND 



We consider a two-component elongated (along the x-direction) repulsive BEC, composed of two different hyperfine states of 
the same alkali isotope. In the case of a highly anisotropic trap (i.e., if the longitudinal and transverse trapping frequencies are 
such that uj x <C oj±), this system can be described by two coupled Gross-Pitaevskii equations (GPEs) of the form d): 

( h 2 2 \ 

Here, ^j(x, t) (j = 1, 2) denote the mean-field wave functions of the two components (normalized to the numbers of atoms 
Nj = \iljj\ 2 dx), m is the atomic mass, and fij are the chemical potentials; furthermore, gjk = 2hu±ajk are the effective 
one-dimensional (ID) coupling constants, denote the three s-wave scattering lengths (note that a\2 — &21) which account 
for collisions between atoms belonging to the same (a,jj) or different (a^, j ^ k) species, while V(x) = (l/2)muj 2 x 2 is the 
external trapping potential. 

Let us now assume that the two-component BEC under consideration consists of two different hyperfine states of 87 Rb, such 
as the states |1, — 1) and |2, 1) used in the experiment of Ref. |30], or the states |1, — 1) and |2, — 2) used in the experiments of 
Refs. QiO ED El. In the first case the scattering lengths take the values an = 100. 4ao, a\2 = 97.66ao and CL22 = 95.00ao, 
while in the second case the respective values are an = 100. 4ao, a\2 = 98.98ao and CL22 = 98.98ao (where ao is the Bohr 
radius). In either case, it is clear that the scattering lengths take approximately the same values, say a^ « a. This way, measuring 
the densities | 2 , length, time and energy in units of 2a, a± = ^/H/uj±, ujJ 1 and Huj±, respectively, we may cast Eqs. M into 
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the following dimensionless form, 

idtut = - X -dlu x +V(x) Ul + (\ Ul \ 2 + \u 2 \ 2 (2) 

id t u 2 = ~d 2 x u 2 + V(x)u 2 + (\ Ul \ 2 + \u 2 \ 2 - n)u2, (3) 

where we have also assumed that the chemical potentials characterizing each component are equal. Note that the potential in 
Eqs. is now given by V(x) — (l/2)Q 2 x 2 , where ft = uj x /uo±_ is a natural small parameter of the system. 

The above system of Eqs. ([2])-([3]) is invariant under SU(2) rotations l28l . In particular, let us first recall that a general matrix 
element of SU(2) takes the form 

L ~ 1 a* 

where a and f3 are complex constants such that \a\ 2 + \/3\ 2 = 1. Then, it can be shown that if {u\,U2) T are solutions of 
Eqs. then, 



U 



U\ \ _ ( OLU\ — fi*U2 



are also solutions of Eqs. Q-Q. This suggests that we may start from the exact dark-bright (DB) soliton solution (which exists 
in the absence of the potential) and, derive the beating dark-dark (DD) soliton solution. More specifically, in the absence of the 
external potential (V(x) = 0), and for the boundary conditions \ui\ 2 — > /i and \u2\ 2 —> as \x\ —> oo, Eqs. ([2])-([3]) possess an 
exact analytical single DB soliton solution of the following form: 

ui(x,i) = ^{cos (j) ta,nh {; + i sin (j)} , (4) 
U2(x,t) = f]sech( i ex.p{ikx + i0(t)}, (5) 

where ^ = D(x — xo(t)), § is the dark soliton's phase angle, cos (j) and r\ represent the amplitude of the dark and bright solitons, 
and D and xq (t) are associated with the inverse width and the center position of the DB soliton. Furthermore, k = D tan <p and 
0(t) are the (constant) wavenumber and phase of the bright soliton, respectively. The above parameters of the DB soliton are 
connected through the following equations: 

D 2 = /icos 2 ^-^ 2 , (6) 
xq = k = D tan 0, (7) 

6 = \(D 2 -k 2 ), (8) 

with xo and denoting the DB soliton velocity and angular frequency, respectively (overdots denote time derivatives). Thus, 
the DB soliton ([?]), §5§ is characterized by three free parameters (seven parameters /i, 0, rj, k, D, xq, 6 and four constraints 
([6])-([8])). Notice that the amplitude r] of the bright soliton, the chemical potential /i of the dark soliton, as well as the (inverse) 
width parameter D of the DB soliton are connected to the number of atoms Nb of the bright soliton by means of the following 
equation: 

N B = f \ U2 \ 2 dx= 2 -y^. (9) 



D 



According to the above arguments, one may start from the DB soliton and construct SU(2) rotated solutions, in the following 
form: 



u\ (x,t) = ay7^{ cos ^ tanh £ + isincj)} — /3*?7 sech^ exp{z/cx + i0(i)}, (10) 
U2(x,i) = /3y7^{ cos ^^anh^ + isin^} + a*^sech^exp{z/cx + i(9(t)}. (11) 

With the additional four parameters a, j3 G C and the constraint \a\ 2 + \ f3\ 2 = 1, the solution ([T0|)-([TT]) is characterized by six 
free parameters. Introducing a new parameter c, the velocity of the background fluid, another solution can be constructed from 
Eqs. ([TQ|)-([TT]) via a Galilean boost: ex.p[i(cx — c 2 t/2)]ui^(x — ct, t). Thus, in the most general case, this DD soliton solution 
is characterized by seven free parameters. One natural set of parameters can be found from the far-field, \x\ —> oo behavior 
consisting of two densities, an overall fluid velocity, and four phases. 
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Due to Galilean invariance and phase invariance, Uj(x,i) = e llfi iUj(x,t), we will assume, without loss of generality, that 
the background is at rest (c = 0) and focus, more specifically, on the case of the SO(2) rotated DB soliton. In this case, the 
corresponding orthogonal matrix is given by: 

u= (cos(x) (12) 



sin(x) cos(x) 

where x is an arbitrary angle. This way, the relevant SO(2) rotated soliton solution takes the form: 

u\(x,t) = cos(x)^/i^{cos(j)t^Lnh(D(x — x o(t))) + is'mcj)} — s'm(x)v sec h(D(x — xo(t))) ex.p{ikx + i6(t)}, (13) 
U2{x,t) = sm(x)^/^^{cos(j)ta 1 nh(D(x — x o(t))) + i s'mcj)} + cos(x)7?secri(Z)(x — xo(t))) ex.p{ikx + i6(t)}, (14) 

The solution (l^-(|14|) is a DD soliton solution characterized by 4 free parameters. The asymptotics of this solution are \ui | 2 — » 
/icos 2 (x) and~]^2|^-> Msin 2 (x) as \x\ — » oo. The densities of the above dark solitons read: 

m = \ui\ 2 = /icos 2 (x) — (/icos 2 (x) cos 2 cj) — rf sin 2 (x))sech 2 £ 

— y7^^ sm (2x) {sin0sin[/cx + ^(t)] + cos^>cos[/cx + ^(t)] tanh^} sech^, (15) 

ri2 = \u2\ 2 = /i sin 2 (x) — (/i sin 2 (x) cos 2 ^) — 7] 2 cos 2 (x))sech 2 ^ 

+ y/jir) sin(2x) {sin sin [/ex + + cos </>cos[/cx + tanh^} sech^, (16) 

while the total density n to t of the DD soliton is given by: 

ntot = ni + n 2 = /i - Z} 2 sech 2 £. (17) 

Notice that the total density of the DD soliton is time-independent and has the form of a dark soliton density of depth D 2 on top 
of a background density /i. The above density is, in fact, identical to the density of the DB soliton; this is due to the fact that 
under SO(2) rotation the total density, as well as all other conserved quantities of the system, remain unchanged. This will be 
particularly important when considering the motion of the DD soliton in a trap — see below. 

On the other hand, one may consider the individual dark soliton densities, n\ and ri2, across the trajectory of the DD soliton, 
i.e., for £ = 0: in such a case, x = xo(t) = kt and the densities read: 



rci(f = 0) = /icos 2 (x)sin 2 ^ + ^ 2 sin 2 (x) 
— \JJir} sin(2x) sin </>sin 



2 

^2(£ = 0) = \± sin 2 (x) sin 2 (j) + rj 2 cos 2 (x) 
+ y/Jirj sin(2x) sin sin 



l (k 2 + D 2 )t 



l(k 2 +D 2 )t 



(18) 



(19) 



It is clear that ni^iS, = 0) are periodic functions of time; the relevant angular frequency (which constitutes the internal beating 
frequency of the DD soliton) is given by: 

^0 = \(k 2 + D 2 ) = - r) 2 sec 2 </>), (20) 

where we have also used Eq. ([6]). The frequency ujq is bounded by two limiting values. First, in the case 77 — » 0, the beating 
DD soliton becomes a regular DD soliton, characterized by a width D = ^JJi cos (j) and a velocity k = yf]lsm.(j)\ in this 
case, ujo — >> (l/2)/i. Second, in the limiting case D — >> 0, the beating DD soliton is reduced to a plane wave; in this case, 
cjo — >■ (1/2) /c 2 . In other words, the intrinsic oscillation frequency take values in the range: 

1 9 1 

-k 2 <UJ < -jJL. (21) 



III. DARK-DARK SOLITONS AS PERIODIC ORBITS IN THE MANAKOV MODEL 



In this section, we analyze the existence, stability and dynamics of single beating DD solitons in a trap of the form V(x) = 
^Q 2 x 2 , considering them as periodic orbits. In the presence of the trap, the dynamics of the center of mass xo(t) of the beating 
DD soliton is still described by the dynamics of the original (unrotated) DB soliton center xo . This is due to the fact that the 
GPEs |2])-([3} are invariant under SO(2) rotations even in the presence of V(x), and so are all conserved quantities of the system, 
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FIG. 2: Left panel: Profiles and densities of the beating dark-dark soliton solution with Q, = 0, ujo = 0.5 and /n = 1.5 at t = 0. Right panel: 
Floquet multipliers spectrum for the dark-dark soliton displayed in the left panel. 



such as the total energy. Since the derivation of the equation of motion for the DB soliton center xq in Ref. |[T6ll was relying on 
the change of energy (due to the presence of the trap), it is clear that the evolution of the beating DD soliton center follows the 
same dynamics: it performs a harmonic oscillation in the trap according to the equation £ + w% sc xo — 0> where the oscillation 
frequency uj osc is given by lfT6l : 

where r = ^= is a measure of the ratio of the number of atoms in the bright and dark soliton component. In order to compute 
the soliton profile and determine its stability, we consider the solution of Eqs. Q-([3]), with gn = #22 = #12 = 1, as a Fourier 
series expansion of period uo l32lL namely, 

CO CO 

uiOM)= E z k(x)e ikulot , u 2 (x,t)= £ y k {x)e ik ^\ (23) 

k= — co k— — oo 

with {z k }, {yk} £ Then, the dynamical equations are reduced to a set of coupled equations: 

\f> ~ k ^0 ~ V(x)]z k + ^plz k = ^2( Z P Z Q + Vpyl) Z k-p+q (24) 

p q 

[/i - ku - V(x)]y k + ^ply k = ^2( z p z q + VpVDVk-p+q (25) 

p q 

where we have used the notation z k = z k (x), y k = y k { x )- If me tra P is absent, it is straightforward to see that 



z o( x ) = \ o tanh(v / 2c^^) = 1/oW, (26) 




^i(x) = -W^ -cj sech( v / 2a^x) = -?/i(x), (27) 



2 

^■(a;) = %(x) = 0, |j|>1 or j = -l, (28) 



is actually the solution ( 13 )-( 14) for x — tt/4, </> = fc = 0, and ujq = D 2 /2. In order to numerically find a DD soliton solution 



in the system with the trap, the previous solution (with the dark component {z k } multiplied by the Thomas-Fermi cloud) is 



introduced as a seed for a fixed-point method in the system of Eqs. ([24]) -([25]). Throughout this section, we have considered — for 
convenience — a trap strength ft = 0.2 in order to consume less time in the numerical calculations. Figures [2] and [3] show the 
periodic orbit for t = without and with a trap potential, respectively. It is worth remarking that solutions in the trap exist for 
/i > 2uo, as predicted in the end of section II. 

Once a periodic solution is found, its (linear) orbital stability can be analyzed by means of Floquet analysis. To this end, 
the time evolution of a small perturbation t), ^(^, t) to a periodic solution {u\${x, t), ^2,0(^5 1)} must be traced. These 
perturbations are introduced in the dynamical equations as: 

ui(x,t) = [ui t o(x,t) + 5£i(x,t)] , u 2 (x,t) = [1/2,0 OM) + *60M)] ^ (29) 
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FIG. 3: Same as Fig. [2]but in the trapped case with Q, = 0.2. 



and the resulting equation at order 0(S) reads: 

1 



idtt.2 



[-^dl + 2|^i 5 o| 2 + |^2,o| 2 - /i + V(x)}^ + uj^l + 1x3,0^1,06 + ^2,0^1,0^2 » 

|^l 5 o| 2 + 2|ii 2 ,o| 2 - /i + V"(x)]^2 + ^2,0^2 + ^1,0^2,06 + ^1,0^2,0^- 



1 



(30) 
(31) 



Then, in the framework of Floquet analysis, the stability properties of periodic orbits are resolved by diagonalizing the mon- 
odromy matrix M, which is defined as: 



Re(£i(z,T)) 
Im(£i(x,T)) 
Re(6(x,T)) 





rRe(^(x,0))-| 


= M 


Hi(i,o)) 


Re(6(x,0)) 

LM6(i,o))J 



(32) 



with T = 2tt/ujo. As the system is symplectic and Hamiltonian, the linear stability of the solutions requires that the monodromy 
eigenvalues, A (also called Floquet multipliers) must lie on the unit circle (see, e.g., (35][36]l for details). The Floquet multipliers 
can also be written as A = exp(zO), with 6 being denoted as the Floquet exponent. An internal mode of the soliton corresponds 
to a spatially localized solution of Eqs. ( [30] )- ( |3T] ), with its oscillation frequency related to the Floquet exponents as uj m — 
Bcjo/(27t). Figures [2] and [3] show a typical Floquet multiplier spectra, indicating stability of the periodic orbits. All the analyzed 
solutions (i.e. with ft — and Q = 0.2) are stable. 

Some interesting results can be extracted by the analysis of the internal modes of the periodic orbits. Figure |4jleft) shows 
the dependence of three internal modes of the Floquet spectrum with respect to /i for ujq = 0.5. The blue line is close to the 



frequency predicted by Eq. ( [22] ) [depicted as a dashed red line] . Indeed, perturbing the beating DD soliton with the corresponding 
eigenmode, we have confirmed that this perturbation leads to an oscillation of the soliton in the trap with a frequency equal to 
that of the eigenmode (cf. left panel of Fig. [5]). It can be observed that the agreement between the numerical eigenfrequency 
and that predicted by Eq. ( [22] ) improves when fi increases, as expected. The right panel of Figure [4] shows the dependence of the 
frequency of the internal mode corresponding to the oscillation of the trap with respect to cjq for fixed /i = 5 and compares it 



with the frequency predicted by Eq. ( [22] ). 

We note here, as an aside in the case Vt = 0, that the internal soliton modes are neutral modes located at (1,0) on the unit 
circle. In particular, the mode associated with the oscillation of the DD soliton in the trap becomes in this case a neutral mode 
associated with the translation of the soliton. The algebraic multiplicity of the multiplier at (1, 0) in the case of Q = is 8, while 
in the trapped case (due to the lifting of translational invariance) it is 6. 

In order to observe the properties of other internal modes, we have perturbed the beating DD soliton with the corresponding 
eigenmodes. In particular, a perturbation along the direction of the localized mode depicted in black in the left panel of Fig. [4] 
leads to a breathing in the width of the soliton — see left panel of Fig. [5] On the other hand, a perturbation along the direction of 
the mode depicted in green in the left panel of Fig. [4] leads to an oscillation of the condensate along the trap, leaving the beating 
DD soliton unaffected (i.e., the soliton stays at the trap center) — see right panel of Fig. [5] 

Finally, we make a remark about the way we have calculated the value Nb that must be introduced in Eq. ( [22] ). The procedure 
consists in performing an SO(2) rotation with x — — tt/4 to the periodic DD soliton at t = 0. This solution is shown in the left 
panel of Fig. [6] whereas the rotated solution is depicted in the right panel of the same figure. Thus, Nb is the norm of the bright 
component of the rotated solution. It can also be inferred from the Fourier coefficients of the periodic orbit: 



N E 



I \u 2 \ 2 dx=\ I [|^^| 2 + |^^| 2 -2Re((^4)(E^))] d ^ 



(33) 
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FIG. 4: (Left panel) Dependence of the eigenfrequencies of some internal modes (solid lines) with respect to \i with fixed cjo = 0.5 (The black, 
blue, green and red lines are the upper solid, middle solid, lower solid and the dashed respectively). The right panel shows that dependence 
with respect to ujq with fixed \i — 5 for the mode in blue in the left panel. In both panels, the red dashed line corresponds to the oscillation 
frequency p2| in a trap with Vt — 0.2. 
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FIG. 5: Contour plots showing the evolution of the densities of the DD soliton components when perturbed by three different eigenmodes: the 
left panel corresponds to a soliton perturbed along the blue mode and leads to trap oscillations; in the central panel, the perturbation (along the 
black mode) leads to a breathing of the soliton width, whereas in the right panel (perturbation along the green mode), the outcome corresponds 
to an oscillation of the whole condensate. In both cases, /x = 3, cjq = 0.5 and Vt — 0.2. 



IV. SINGLE BEATING DARK-DARK SOLITONS NEAR AND FAR FROM THE MANAKOV LIMIT 



We now turn to a numerical study of the properties of the beating DD soliton states. Firstly, in the absence of a trap, we 
are going to compare the integrable case with equal scattering lengths gn : gi2 : #22 = 1 • 1 ' 1 to the non-integrable case 
9u • 9 12 • 922 = 1.03 : 1 : 0.97 (see Ref. [ 30 1). From Fig. [7] we observe that both of the dark components are oscillating with 
fixed frequencies and these two cases are very similar l38l (a feature that we have generically observed between the dynamical 
phenomenology of these two cases). 

To highlight the fact that substantial variations of the scattering length-which can be imposed by virtue of a Feshbach 
resonance-may have a significant impact on the robustness of the beating DD solitons, we consider scattering lengths in the 
set with ratios gn : #12 : #22 = 9 ' 1 : 1- In particular, we take g = 1.1, 1.2, 1.6 in Fig. [8] When g is not so large, i.e., 
g = 1.1, 1.2, the beating DD soliton oscillates and, as t increases, the change in the scattering length results in mobility of 
the coherent structure. However, more dramatic events can arise when g is relatively large, e.g., for g = 1.6. There, we can 
see that the soliton is finally split into two fragments (upon growth of the intrinsic beating oscillation which eventually induces 
the splitting) and results in two states that resemble dark-antidark solitons [TT] (see also Ref. |29l ). In particular, each of the 
components acquires a dark soliton coupled to a lump in the second component. 

In Fig. [9] we show a particular example of the DD soliton in the trap, which oscillates around its center; the parameter values 
are /i = 1, 77 = 0.6, initial soliton position xo(t = 0) = 2.5, and trap strength Q = 0.05. Note that for these runs, the initial 
profile of the beating DD soliton in the trap is approximated by the numerically found (in trap) ground state — i.e., the Thomas- 
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FIG. 6: (Left panel) Profiles and densities of a periodic orbit at t = with \i — 3, cjo = 0.5 and Q = 0.2. The right panel shows the 
dark-bright soliton arising by rotating with \ — — tt/4 the dark-dark soliton of the left panel. 





FIG. 7: The comparison between the integrable case gn : gi2 : #22 = 1:1:1 (left column) and the non-integrable case gn : gi2 : #22 = 
1.03 : 1 : 0.97 (right column), is demonstrated. The upper panels show the densities of the first dark component while the lower ones show 
the second dark component. Here 77 = 0.6, x — tt/4, = 0, k = 0, \i — 1. Based on the similarity of the relevant dynamics, we will focus on 
the case of unit nonlinear coefficients. 



Fermi cloud — multiplied by the beating DD solution (without a trap) of Eqs. ([T3])-([T4|). Then via a time-stepping algorithm 
(a fourth-order Runge-Kutta scheme), we obtain the time evolution of the densities of the oscillating solitons in the upper two 
panels. Moreover, the left lower panel shows the center of mass of the beating DD soliton in the trap. Using Fourier analysis, 
we can infer the numerical frequency of in-trap oscillation, which can, in turn, be compared to the analytical one, cf. Eq. (22). 
As shown in the bottom right panel of the figure, there is very good agreement between the two. 

Next, we consider the in-trap dynamics of a single beating DD soliton but for the non-integrable cases. Again, when gn : 
9 12 922 = 1-03 : 1 : 0.97, we observe a nearly identical phenomenology to that of unity's. For the more significant deviations 
from that case of the form gn : gi2 : #22 = 9 ' 1 : 1 where g = 1.1, 1.2, 1.6, the results are reported in Fig. 10 For lower values 
of g = 1.1, 1.2, the behavior of the DD is similar to the case with g = 1, however, we progressively observe more significant 
radiative emissions which also affect the oscillation frequency. However, once again the modifications of the phenomenology 
are most dramatic in the case of g = 1.6 of the bottom panels. There, the radiation emission is accompanied by growing intrinsic 
oscillations which eventually result in the breakup and formation of a single dark-antidark solitary wave. 



V. TWO DARK-DARK SOLITON STATES: DYNAMICS AND INTERACTIONS 

We now consider the interactions of two beating DD solitons. We once again start from the untrapped case and use as an 
initial ansatz a two-DB soliton state of the form: 

u\ = (cos </>tanh£_ + i sin <p) (cos </>tanh£ + — i sin <j)) (34) 
u 2 = f]^H-e i(kx+m) + e iA %secK + e i( - fc ^ (£)) (35) 
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FIG. 8: The dynamics of the densities of a DD soliton in the absence of a trap but for # = 1.1, g = 1.2, g = 1.6 respectively (g = #n). When 
g = 1.1 or g = 1.2, the soliton is set into translational motion. However, for g = 1.6 the coherent structure executes a growing oscillation 
which eventually results in its splitting into a pair of dark-antidark solitons (i.e., a dark soliton in the one component coupled to a lump in the 
other). The parameters used here are the same as in Fig. [7] 
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FIG. 9: Motion of a dark-dark soliton in a trap of strength Q = 0.05. Other parameters are the same as in Fig. [7] The upper two panels 
show the oscillation of the wave initially centered at xo = 2.5 (the chemical potential is /i = 1). The lower left panel demonstrates the center 
of mass of the DD in the upper panels. The analytical oscillation frequency, given by Eq. {22} , is 0.03123, while the numerical frequency, 
calculated by Fourier transform, is 0.03238. The lower right panel yields the comparison between the analytically calculated frequencies (red 
line) versus the numerical obtained ones (the blue triangles), as r varies between 0.1 and 14. 



where £± = D{x ± xq), 2xo is the relative distance between the two solitons, and AO is the relative phase between the two 
bright solitons. Below we consider both the out-of-phase (OOP) case, AO — tt, as well as the in-phase (IP) case AO = 0. Once 
again taking advantage of the model invariance under the SO(2) rotations, as we did for the single DD soliton case, we use the 
orthogonal matrix $V2) and obtain a two-beating-DD- soliton ansatz in the form: 

ui = cos(x) (cos (/>tanh£_ + zsin0) (cos (/>tanh£ + — ismcf)) 

- sin(x) (r ] sech£_e i ( kx + m) + e iA ^secK + e i( -^ + ^ (t)) ) , (36) 

U2 = sin(x) (cos </> tanh + isin0) (cos0tanh£ + — ismc/)) 

+ cos(x) (77sech£_e^+^ + e iAG 7]sech^- kx+G ^ . (37) 



In our numerical study for the dynamics of the two-beating-DD- soliton state, we first consider the integrable case, correspond- 
ing to = #12 = #22 = 1, both for the in-phase and out-of-phase cases. The results of the simulations, using initial conditions 
corresponding to the above ansatz, are shown in Fig.[TT] In the in-phase case, the repulsion between the beating DD solitons is 
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FIG. 10: The comparison of the oscillation of the density of a DD soliton within a trap of trap frequency Q — 0.05 for different values of 
g. The soliton is initialized at xq — 2.5; g is set to be 1.1 (top panels), 1.2 (middle panels), 1.6 (bottom panels) for the combination of the 
scattering length gu \ g\2 '■ #22 = g : 1 : 1. Other parameters are similar to Fig. [9] 
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FIG. 11: Space-time contour plots of two beating DD soliton densities in phase (left) and out of phase (right) for gu — gi2 — #22 = 1 . Here 

x = tt/4, 77 = 0.5, xo = 1.5, D = 1.2. 



immediately evident resulting in the strong separation of the two waves (which still perform their internal beating). On the other 
hand, in the out-of-phase case, the competition between the repulsion of the dark components and the attraction between the 
bright components of the progenitor DB solitons (see Ref. l33l ) can be discerned, as the configuration remains nearly stationary 
for a lengthy evolution interval. Finally, however, the repulsive interaction prevails and the solitons eventually separate. 

Next, we consider the non-integrable case. Since for gu : g 12 : #22 = 1-03 : 1 : 0.97, the phenomenology is again very 
similar to gu = g\2 — #22, we consider the significant departure from this limit pertaining to gu : gi2 : #22 = 1.6 : 1 : 1. 
In Fig. [T2j we observe that in the in-phase case, the two beating DD solitons initially separate and move away from each other, 
then they are reflected from the domain boundary and a new collision occurs. After this collision, a highly nontrivial event is 
observed, namely one of the two beating DD solitons is decomposed into a dark-antidark soliton pair, with each of these solitons 
moving with different velocities. For the out-of-phase case, the separation arises much faster than for the unit coefficients and, 
interestingly, results in an asymmetric evolution with one of the DD solitons breaking up in a pair of dark-antidark solitons (as 
in Fig.[S]of section II). Notice that, as in the in-phase case, the other soliton is not broken up in a similar way during the horizon 
of the simulation although it is likely that such an event will also occur for that wave. 

Next, we consider the two-beating-DD soliton in the trap, in the case of unit coefficients. We set V(x) — \Vt 2 x 2 , with 



Q = 0.05, and the chemical potential [i = 1. From Fig. [T3J we infer that the two beating DD solitons are now trapped and 
oscillate around an equilibrium position. Notice that in the in-phase case, the solitons perform out-of-phase oscillations and 
undergo quasi-elastic collisions, while In the out-of-phase case, the weak residual repulsion is counter balanced by the presence 
of the trap, and we observe that the two beating DD solitons remain in a close distance to each other. 
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FIG. 12: Space-time contour plots of two DD soliton densities in phase (left) and out of phase (right) for g = 1.6 in the set gn : g±2 : #22 
1.6 : 1 : 1. Here x = tt/4, 77 = 0.5, x = 1.5, D = 1.2. 
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FIG. 13: Space-time contour plots of two beating DD soliton densities in phase (left) and out of phase (right) in the case of equal gij's within 
a harmonic trap with trap frequency Q = 0.05. Here x — ?r/4, rj — 0.5, xq = 1.5, D = 1.2, /x = 1. 



Finally, we consider two DD with gu : g\2 ' #22 = 1.6 : 1 : 1 within the same trap. In this case, we observe that despite 
the presence of the trap, it is not possible to sustain a robust set of oscillations and interactions between the beating DD solitons. 
This is especially true in the out-of-phase case, where the oscillatory growth of the beating eventually results in the breakup of 
the DD soliton states into dark-antidark ones (which generally appear more robust for such higher values of g). 



VI. CONCLUSIONS AND FUTURE CHALLENGES 



In this work, we have studied the stability and dynamics of beating dark-dark (DD) solitons in pseudo-spinor Bose-Einstein 
condensates, motivated by recent experiments where such structures were observed. We have illustrated the connection of these 
solitons with internal density oscillations to dark-bright (DB) solitons identified earlier, through SO(2) (and more generally 
SU(2)) rotations. We have illustrated that such states persist in the presence of the trap and, in fact, oscillate with the frequency 
previously predicted for dark-bright solitons. Using Floquet analysis, we have also identified beating dark-dark solitons as stable 
periodic orbits in the integrable (Manakov) limit with and without a trap. 

We have also investigated in detail the effect of the deviation from the Manakov case by considering different from unity 
scattering length ratios. We have shown that when the deviation from the integrable case is small (as is the physically relevant 
case of a pseudo-spinor condensate composed by different spin states of rubidium), then the stability and dynamics of beating 
dark-dark solitons follow that of the integrable case. However, we also illustrated that a significant departure of the ratios 
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FIG. 14: Space-time contour plots of two beating DD soliton densities in phase (left) and out of phase (right) for the case with gn : gi 2 : 
#22 = 1.6 : 1 : 1 within a harmonic trap with trap frequency Q = 0.05. The parameters used are the same as for the previous figure. 



of the scattering lengths from this limit (towards the miscible regime) will eventually break up beating dark-dark solitons in 
favor of dark-antidark soliton entities. We have also considered the interaction of beating dark-dark solitons finding a typically 
repulsive dynamical behavior, which can be attenuated only in the case where the bright components (of the progenitor dark- 
bright solitons, used to create the dark-dark ones) are out-of-phase (and, hence, attracting each other). In that case, especially in 
the presence of a trap, a robust set of multi-beating-dark-dark- soliton states can be created. 

The discussion of DD solitons in this work has focused upon those states that can be constructed, in the spatially extended, 
Manakov case, from the SU(2) rotation of a DB soliton and confined states in the presence of a trapping potential. In both cases, 
each component of the DD soliton exhibits the same background flow velocity. In a series of experiments (l8l[T9l[25j[33l, a 
relative flow between two condensate components induced by a magnetic field gradient led to DB solitons and counterflow- 
induced modulational instability resulting in the formation of a number of beating DD solitons. It is natural, then, to inquire 
into the effect that relative motion between two condensate components has on localized structures. In the integrable case, the 
most general DD soliton was constructed using a Backhand transformation [28]. Because it allows for a counterflow, this soliton 
is characterized by eight free parameters in contrast to the seven parameter SU(2) rotated DB soliton studied here or the seven 
parameter static DD soliton [37]. However, and since the present study focused predominantly on the five-parameter family 
stemming from the SO(2) rotation, the persistence, stability, and interactions of the full seven parameter solitonic states (and 
even the eight parameter generalization thereof presented within [28]) in the non-integrable case constitute themes worthy of 
further study. 

There are many other directions that are worth considering further along the lines of this work. Quantifying further (and 
semi-analytically, if possible) the interactions between the beating dark-dark solitons, as well as studying in more detail the 
dark-antidark solitons that appear to spontaneously arise from their breakup in the miscible regime are interesting extensions of 
this work in the one-dimensional setting. On the other hand, one naturally may consider the two-dimensional (2D) generalization 
of the considerations herein, especially upon bearing in mind that the SU(2) (or SO(2)) rotations used herein are not restricted to 
the one-dimensional realm in any particular way. In that regard, one may envision vortex-bright soliton states l34l (i.e., the 2D 
analog of the dark-bright waves) rotated via SO(2) to produce vortex- vortex type states (in analogy to the dark-dark ones). Such 
states are currently under study and will be reported in future publications. 
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